These slides and all other course materials can be found at
github.com/brry/rhydro #slides
To get all the material including the datasets and presentation source code, we recommend to download the whole github course repository.

This is an R Markdown Notebook.
For discussions, please visit the Hydrology in R Facebook group.
Before running the code blocks below, we suggest to get package installation instructions by running:

source("https://raw.githubusercontent.com/brry/rhydro/master/checkpc.R")


Aim and contents of this workshop

We want to:

We can not:

We have prepared:

 

Before we get started, please let us know your current R knowledge level by filling out the short survey at
bit.ly/knowR


top

Report

Good coding practice, report generation (Rstudio, rmarkdown, R notebook)
Daniel Klotz


top

GIS

Using R as GIS (reading a rainfall shapefile + Kriging, sf + leaflet + mapview + OSMscale)
Berry Boessenkool

Shapefiles

Reading shapefiles with maptools::readShapeSpatial and rgdal::readOGR is obsolete.
Instead, use sf::st_read. sf is on CRAN since oct 2016.
Main reaction when using sf: “Wow, that is fast!”
Download the shapefile or better: download the whole github course repository

rain <- sf::st_read("data/PrecBrandenburg/niederschlag.shp")
Reading layer `niederschlag' from data source `/home/berry/Dropbox/R/rhydro/presentations/data/PrecBrandenburg/niederschlag.shp' using driver `ESRI Shapefile'
converted into: POLYGON
Simple feature collection with 277 features and 1 field
geometry type:  POLYGON
dimension:      XY
bbox:           xmin: 3250175 ymin: 5690642 xmax: 3483631 ymax: 5932731
epsg (SRID):    NA
proj4string:    +proj=tmerc +lat_0=0 +lon_0=15 +k=0.9996 +x_0=3500000 +y_0=0 +ellps=GRS80 +units=m +no_defs
centroids <- sf::st_centroid(rain)
centroids <- sf::st_coordinates(centroids)

top

Plotting, maps

Static plot:

plot(rain[,1])

Static map:

prj <- sf::st_crs(rain)$proj4string
centroids <- as.data.frame(centroids)
cent_ll <- OSMscale::projectPoints(Y,X, data=centroids, to=OSMscale::pll(), from=prj)
map_static <- OSMscale::pointsMap(y,x, cent_ll, fx=0.08, type="maptoolkit-topo", zoom=6)
Downloading map with extend 11.029332, 14.981464, 51.10418, 53.765423 ...
Done. Now plotting...

Interactive map:

library(leaflet)
cent_ll$info <- paste0(sample(letters,nrow(cent_ll),TRUE), ", ", round(cent_ll$x,2), 
                       ", ", round(cent_ll$y,2))
leaflet(cent_ll) %>% addTiles() %>% addCircleMarkers(lng=~x, lat=~y, popup=~info)

Interactive map of shapefile:

# devtools::install_github("environmentalinformatics-marburg/mapview", ref = "develop")
library(berryFunctions) # classify, seqPal
col <- seqPal(n=100, colors=c("red","yellow","blue"))[classify(rain$P1)$index]
mapview::mapview(rain, col.regions=col)

top

Kriging

Plot original points colored by third dimension:

pcol <- colorRampPalette(c("red","yellow","blue"))(50)
x <- centroids$X # use cent_ll$x for projected data
y <- centroids$Y
berryFunctions::colPoints(x, y, rain$P1, add=FALSE, col=pcol)

Calculate the variogram and fit a semivariance curve

library(geoR)
--------------------------------------------------------------
 Analysis of Geostatistical Data
 For an Introduction to geoR go to http://www.leg.ufpr.br/geoR
 geoR version 1.7-5.2 (built on 2016-05-02) is now loaded
--------------------------------------------------------------
geoprec <- as.geodata(cbind(x,y,rain$P1))
vario <- variog(geoprec, max.dist=130000) # other maxdist for lat-lon data
variog: computing omnidirectional variogram
fit <- variofit(vario)
variofit: covariance model used is matern 
variofit: weights used: npairs 
variofit: minimisation function used: optim 
initial values not provided - running the default search
variofit: searching for best initial value ... selected values:
              sigmasq   phi        tausq kappa
initial.value "1325.81" "19999.05" "0"   "0.5"
status        "est"     "est"      "est" "fix"
loss value: 104819060.429841 
plot(vario)
lines(fit)

Determine a useful resolution (keep in mind that computing time rises exponentially with grid size)

# distance to closest other point:
d <- sapply(1:length(x), function(i)
            min(berryFunctions::distance(x[i], y[i], x[-i], y[-i])) )
# for lat-long data use (2017-Apr only available in github version of OSMscale)
# d <- OSMscale::maxEarthDist(y,x, data=cent_ll, fun=min)
hist(d/1000, breaks=20, main="distance to closest gauge [km]")

mean(d/1000) # 8 km
[1] 8.165713

Perform kriging on a grid with that resolution

res <- 1000 # 1 km, since stations are 8 km apart on average
grid <- expand.grid(seq(min(x),max(x),res),
                    seq(min(y),max(y),res))
krico <- krige.control(type.krige="OK", obj.model=fit)
#krobj <- krige.conv(geoprec, loc=grid, krige=krico)
#save(krobj, file="data/krobj.Rdata")
load("data/krobj.Rdata") # line above is too slow for recreation each time

Plot the interpolated values with or an equivalent (see Rclick 4.15) and add contour lines.

par(mar=c(0,3,0,3))
geoR:::image.kriging(krobj, col=pcol)
colPoints(x, y, rain$P1, col=pcol, legargs=list(horiz=F, title="Prec",y1=0.1,x1=0.9))
points(x,y)
plot(rain, col=NA, add=TRUE)


top

Discharge

River discharge time-series visualisation and extreme value statistics (animation + extremeStat)
Berry Boessenkool


top

Hydmod

Hydrological modelling with airGR
Katie Smith


top

EDA

Exploratory Data Analysis including flow duration curve and trend analysis on time-series
Shaun Harrigan


top

Discussion

Please give us feedback at bit.ly/feedbackR

For discussions, please use the Hydrology in R Facebook group.

LS0tCnRpdGxlOiAiVXNpbmcgUiBpbiBoeWRyb2xvZ3kgKEVHVTIwMTcgc2hvcnQgY291cnNlKSIKb3V0cHV0OiAKICBodG1sX25vdGVib29rOiAKICAgIGNvZGVfZm9sZGluZzogbm9uZQogICAgdG9jOiB5ZXMKLS0tCgpUaGVzZSBzbGlkZXMgYW5kIGFsbCBvdGhlciBjb3Vyc2UgbWF0ZXJpYWxzIGNhbiBiZSBmb3VuZCBhdCAgCjxmb250IHNpemU9IjYiPltnaXRodWIuY29tL2Jycnkvcmh5ZHJvXShodHRwczovL2dpdGh1Yi5jb20vYnJyeS9yaHlkcm8pPC9mb250PiAKPGZvbnQgc2l6ZT0iNCI+I3NsaWRlczwvZm9udD4gICAKVG8gZ2V0IGFsbCB0aGUgbWF0ZXJpYWwgaW5jbHVkaW5nIHRoZSBkYXRhc2V0cyBhbmQgcHJlc2VudGF0aW9uIHNvdXJjZSBjb2RlLCB3ZSByZWNvbW1lbmQgdG8KW2Rvd25sb2FkIHRoZSB3aG9sZSBnaXRodWIgY291cnNlIHJlcG9zaXRvcnldKGh0dHBzOi8vZ2l0aHViLmNvbS9icnJ5L3JoeWRyby9hcmNoaXZlL21hc3Rlci56aXApLgoKVGhpcyBpcyBhbiBbUiBNYXJrZG93biBOb3RlYm9va10oaHR0cDovL3JtYXJrZG93bi5yc3R1ZGlvLmNvbS9yX25vdGVib29rcy5odG1sKS4gIApGb3IgZGlzY3Vzc2lvbnMsIHBsZWFzZSB2aXNpdCB0aGUgCltIeWRyb2xvZ3kgaW4gUiBGYWNlYm9vayBncm91cF0oaHR0cHM6Ly93d3cuZmFjZWJvb2suY29tL2dyb3Vwcy8xMTMwMjE0Nzc3MTIzOTA5LykuICAKQmVmb3JlIHJ1bm5pbmcgdGhlIGNvZGUgYmxvY2tzIGJlbG93LCB3ZSBzdWdnZXN0IHRvIGdldCBwYWNrYWdlIGluc3RhbGxhdGlvbiBpbnN0cnVjdGlvbnMgYnkgcnVubmluZzoKYGBgUgpzb3VyY2UoImh0dHBzOi8vcmF3LmdpdGh1YnVzZXJjb250ZW50LmNvbS9icnJ5L3JoeWRyby9tYXN0ZXIvY2hlY2twYy5SIikKYGBgCgpcCgoqKkFpbSBhbmQgY29udGVudHMgb2YgdGhpcyB3b3Jrc2hvcCoqCgpXZSB3YW50IHRvOiAgCgoqIFNob3cgb2ZmIGhvdyBhd2Vzb21lIFIgaXMgZm9yIGh5ZHJvbG9neSAoaXQncyBSLXNvbWUhXl4pICAKKiBDb252aW5jZSB5b3UgdG8gc3RhcnQgb3IgY29udGludWUgdXNpbmcgUiAgCiogUHJvdmlkZSBhbGwgdGhlIGNvZGUgZm9yIHlvdSBhcyBhIHN0YXJ0aW5nIHBvaW50CgpXZSBjYW4gbm90OiAgCgoqIFRlYWNoIHlvdSBhY3R1YWwgUiBjb2RpbmcgKDkwIG1pbnMgaXMgdG9vIHNob3J0IGZvciBhIHR1dG9yaWFsKQoKV2UgaGF2ZSBwcmVwYXJlZDoKCiogW0dvb2QgY29kaW5nIHByYWN0aWNlLCByZXBvcnQgZ2VuZXJhdGlvbl0oI3JlcG9ydCkgKFJzdHVkaW8sIGBybWFya2Rvd25gLCBSIG5vdGVib29rKQoqIFtVc2luZyBSIGFzIEdJU10oI2dpcykgKHJlYWRpbmcgYSByYWluZmFsbCBzaGFwZWZpbGUgKyBLcmlnaW5nLCBgc2ZgICsgYGxlYWZsZXRgICsgYG1hcHZpZXdgICsgYE9TTXNjYWxlYCkKKiBbUml2ZXIgZGlzY2hhcmdlIHRpbWUtc2VyaWVzXSgjZGlzY2hhcmdlKSB2aXN1YWxpc2F0aW9uIGFuZCBleHRyZW1lIHZhbHVlIHN0YXRpc3RpY3MgKGBhbmltYXRpb25gICsgYGV4dHJlbWVTdGF0YCkKKiBbSHlkcm9sb2dpY2FsIG1vZGVsbGluZ10oI2h5ZG1vZCkgd2l0aCBgYWlyR1JgCiogW0V4cGxvcmF0b3J5IERhdGEgQW5hbHlzaXNdKCNlZGEpIGluY2x1ZGluZyBmbG93IGR1cmF0aW9uIGN1cnZlIGFuZCB0cmVuZCBhbmFseXNpcyBvbiB0aW1lLXNlcmllcwoKXCAKCkJlZm9yZSB3ZSBnZXQgc3RhcnRlZCwgcGxlYXNlIGxldCB1cyBrbm93IHlvdXIgY3VycmVudCBSIGtub3dsZWRnZSBsZXZlbCBieSBmaWxsaW5nIG91dCB0aGUgc2hvcnQgc3VydmV5IGF0ICAKPGZvbnQgc2l6ZT0iNiI+W2JpdC5seS9rbm93Ul0oaHR0cHM6Ly9iaXQubHkva25vd1IpPC9mb250PiAKClwKClt0b3BdKCN0b3ApCgojIFJlcG9ydApHb29kIGNvZGluZyBwcmFjdGljZSwgW3JlcG9ydCBnZW5lcmF0aW9uXSgjcmVwb3J0KSAoUnN0dWRpbywgYHJtYXJrZG93bmAsIFIgbm90ZWJvb2spICAKKipEYW5pZWwgS2xvdHoqKgoKXApbdG9wXSgjdG9wKQoKIyBHSVMKVXNpbmcgUiBhcyBHSVMgKHJlYWRpbmcgYSByYWluZmFsbCBzaGFwZWZpbGUgKyBLcmlnaW5nLCBgc2ZgICsgYGxlYWZsZXRgICsgYG1hcHZpZXdgICsgYE9TTXNjYWxlYCkgIAoqKkJlcnJ5IEJvZXNzZW5rb29sKioKCiMjIyBTaGFwZWZpbGVzCgpSZWFkaW5nIHNoYXBlZmlsZXMgd2l0aCBgbWFwdG9vbHM6OnJlYWRTaGFwZVNwYXRpYWxgIGFuZCBgcmdkYWw6OnJlYWRPR1JgIGlzIG9ic29sZXRlLiAgCkluc3RlYWQsIHVzZSBgc2Y6OnN0X3JlYWRgLiBgc2ZgIGlzIG9uIENSQU4gc2luY2Ugb2N0IDIwMTYuICAKTWFpbiByZWFjdGlvbiB3aGVuIHVzaW5nIHNmOiAiV293LCB0aGF0IGlzIGZhc3QhIiAgCltEb3dubG9hZCB0aGUgc2hhcGVmaWxlXShodHRwczovL21pbmhhc2thbWFsLmdpdGh1Yi5pby9Eb3duR2l0LyMvaG9tZT91cmw9aHR0cHM6Ly9naXRodWIuY29tL2Jycnkvcmh5ZHJvL3RyZWUvbWFzdGVyL3ByZXNlbnRhdGlvbnMvZGF0YS9QcmVjQnJhbmRlbmJ1cmcpIApvciBiZXR0ZXI6IFtkb3dubG9hZCB0aGUgd2hvbGUgZ2l0aHViIGNvdXJzZSByZXBvc2l0b3J5XShodHRwczovL2dpdGh1Yi5jb20vYnJyeS9yaHlkcm8vYXJjaGl2ZS9tYXN0ZXIuemlwKQoKYGBge3J9CnJhaW4gPC0gc2Y6OnN0X3JlYWQoImRhdGEvUHJlY0JyYW5kZW5idXJnL25pZWRlcnNjaGxhZy5zaHAiKQpjZW50cm9pZHMgPC0gc2Y6OnN0X2NlbnRyb2lkKHJhaW4pCmNlbnRyb2lkcyA8LSBzZjo6c3RfY29vcmRpbmF0ZXMoY2VudHJvaWRzKQpgYGAKClt0b3BdKCN0b3ApCgojIyMgUGxvdHRpbmcsIG1hcHMKClN0YXRpYyBwbG90OgpgYGB7cn0KcGxvdChyYWluWywxXSkKYGBgCgpTdGF0aWMgbWFwOgpgYGB7cn0KcHJqIDwtIHNmOjpzdF9jcnMocmFpbikkcHJvajRzdHJpbmcKY2VudHJvaWRzIDwtIGFzLmRhdGEuZnJhbWUoY2VudHJvaWRzKQpjZW50X2xsIDwtIE9TTXNjYWxlOjpwcm9qZWN0UG9pbnRzKFksWCwgZGF0YT1jZW50cm9pZHMsIHRvPU9TTXNjYWxlOjpwbGwoKSwgZnJvbT1wcmopCm1hcF9zdGF0aWMgPC0gT1NNc2NhbGU6OnBvaW50c01hcCh5LHgsIGNlbnRfbGwsIGZ4PTAuMDgsIHR5cGU9Im1hcHRvb2xraXQtdG9wbyIsIHpvb209NikKYGBgCgpJbnRlcmFjdGl2ZSBtYXA6CmBgYHtyfQpsaWJyYXJ5KGxlYWZsZXQpCmNlbnRfbGwkaW5mbyA8LSBwYXN0ZTAoc2FtcGxlKGxldHRlcnMsbnJvdyhjZW50X2xsKSxUUlVFKSwgIiwgIiwgcm91bmQoY2VudF9sbCR4LDIpLCAKICAgICAgICAgICAgICAgICAgICAgICAiLCAiLCByb3VuZChjZW50X2xsJHksMikpCmxlYWZsZXQoY2VudF9sbCkgJT4lIGFkZFRpbGVzKCkgJT4lIGFkZENpcmNsZU1hcmtlcnMobG5nPX54LCBsYXQ9fnksIHBvcHVwPX5pbmZvKQpgYGAKCkludGVyYWN0aXZlIG1hcCBvZiBzaGFwZWZpbGU6CmBgYHtyfQojIGRldnRvb2xzOjppbnN0YWxsX2dpdGh1YigiZW52aXJvbm1lbnRhbGluZm9ybWF0aWNzLW1hcmJ1cmcvbWFwdmlldyIsIHJlZiA9ICJkZXZlbG9wIikKbGlicmFyeShiZXJyeUZ1bmN0aW9ucykgIyBjbGFzc2lmeSwgc2VxUGFsCmNvbCA8LSBzZXFQYWwobj0xMDAsIGNvbG9ycz1jKCJyZWQiLCJ5ZWxsb3ciLCJibHVlIikpW2NsYXNzaWZ5KHJhaW4kUDEpJGluZGV4XQptYXB2aWV3OjptYXB2aWV3KHJhaW4sIGNvbC5yZWdpb25zPWNvbCkKYGBgCgpbdG9wXSgjdG9wKQoKIyMjIEtyaWdpbmcKClBsb3Qgb3JpZ2luYWwgcG9pbnRzIGNvbG9yZWQgYnkgdGhpcmQgZGltZW5zaW9uOgpgYGB7cn0KcGNvbCA8LSBjb2xvclJhbXBQYWxldHRlKGMoInJlZCIsInllbGxvdyIsImJsdWUiKSkoNTApCnggPC0gY2VudHJvaWRzJFggIyB1c2UgY2VudF9sbCR4IGZvciBwcm9qZWN0ZWQgZGF0YQp5IDwtIGNlbnRyb2lkcyRZCmJlcnJ5RnVuY3Rpb25zOjpjb2xQb2ludHMoeCwgeSwgcmFpbiRQMSwgYWRkPUZBTFNFLCBjb2w9cGNvbCkKYGBgCgpDYWxjdWxhdGUgdGhlIHZhcmlvZ3JhbSBhbmQgZml0IGEgc2VtaXZhcmlhbmNlIGN1cnZlCmBgYHtyfQpsaWJyYXJ5KGdlb1IpCmdlb3ByZWMgPC0gYXMuZ2VvZGF0YShjYmluZCh4LHkscmFpbiRQMSkpCnZhcmlvIDwtIHZhcmlvZyhnZW9wcmVjLCBtYXguZGlzdD0xMzAwMDApICMgb3RoZXIgbWF4ZGlzdCBmb3IgbGF0LWxvbiBkYXRhCmZpdCA8LSB2YXJpb2ZpdCh2YXJpbykKcGxvdCh2YXJpbykKbGluZXMoZml0KQpgYGAKCkRldGVybWluZSBhIHVzZWZ1bCByZXNvbHV0aW9uIAooa2VlcCBpbiBtaW5kIHRoYXQgY29tcHV0aW5nIHRpbWUgcmlzZXMgZXhwb25lbnRpYWxseSB3aXRoIGdyaWQgc2l6ZSkKYGBge3J9CiMgZGlzdGFuY2UgdG8gY2xvc2VzdCBvdGhlciBwb2ludDoKZCA8LSBzYXBwbHkoMTpsZW5ndGgoeCksIGZ1bmN0aW9uKGkpCiAgICAgICAgICAgIG1pbihiZXJyeUZ1bmN0aW9uczo6ZGlzdGFuY2UoeFtpXSwgeVtpXSwgeFstaV0sIHlbLWldKSkgKQojIGZvciBsYXQtbG9uZyBkYXRhIHVzZSAoMjAxNy1BcHIgb25seSBhdmFpbGFibGUgaW4gZ2l0aHViIHZlcnNpb24gb2YgT1NNc2NhbGUpCiMgZCA8LSBPU01zY2FsZTo6bWF4RWFydGhEaXN0KHkseCwgZGF0YT1jZW50X2xsLCBmdW49bWluKQpoaXN0KGQvMTAwMCwgYnJlYWtzPTIwLCBtYWluPSJkaXN0YW5jZSB0byBjbG9zZXN0IGdhdWdlIFtrbV0iKQptZWFuKGQvMTAwMCkgIyA4IGttCmBgYAoKUGVyZm9ybSBrcmlnaW5nIG9uIGEgZ3JpZCB3aXRoIHRoYXQgcmVzb2x1dGlvbiAKYGBge3J9CnJlcyA8LSAxMDAwICMgMSBrbSwgc2luY2Ugc3RhdGlvbnMgYXJlIDgga20gYXBhcnQgb24gYXZlcmFnZQpncmlkIDwtIGV4cGFuZC5ncmlkKHNlcShtaW4oeCksbWF4KHgpLHJlcyksCiAgICAgICAgICAgICAgICAgICAgc2VxKG1pbih5KSxtYXgoeSkscmVzKSkKa3JpY28gPC0ga3JpZ2UuY29udHJvbCh0eXBlLmtyaWdlPSJPSyIsIG9iai5tb2RlbD1maXQpCiNrcm9iaiA8LSBrcmlnZS5jb252KGdlb3ByZWMsIGxvYz1ncmlkLCBrcmlnZT1rcmljbykKI3NhdmUoa3JvYmosIGZpbGU9ImRhdGEva3JvYmouUmRhdGEiKQpsb2FkKCJkYXRhL2tyb2JqLlJkYXRhIikgIyBsaW5lIGFib3ZlIGlzIHRvbyBzbG93IGZvciByZWNyZWF0aW9uIGVhY2ggdGltZQpgYGAKClBsb3QgdGhlIGludGVycG9sYXRlZCB2YWx1ZXMgd2l0aCBccmNvZGV7aW1hZ2V9IG9yIGFuIGVxdWl2YWxlbnQgCihzZWUgW1JjbGlja10oaHR0cHM6Ly9naXRodWIuY29tL2JycnkvcmNsaWNrKSA0LjE1KSBhbmQgYWRkIGNvbnRvdXIgbGluZXMuCmBgYHtyfQpwYXIobWFyPWMoMCwzLDAsMykpCmdlb1I6OjppbWFnZS5rcmlnaW5nKGtyb2JqLCBjb2w9cGNvbCkKY29sUG9pbnRzKHgsIHksIHJhaW4kUDEsIGNvbD1wY29sLCBsZWdhcmdzPWxpc3QoaG9yaXo9RiwgdGl0bGU9IlByZWMiLHkxPTAuMSx4MT0wLjkpKQpwb2ludHMoeCx5KQpwbG90KHJhaW4sIGNvbD1OQSwgYWRkPVRSVUUpCmBgYAoKXApbdG9wXSgjdG9wKQoKIyBEaXNjaGFyZ2UKUml2ZXIgZGlzY2hhcmdlIHRpbWUtc2VyaWVzIHZpc3VhbGlzYXRpb24gYW5kIGV4dHJlbWUgdmFsdWUgc3RhdGlzdGljcyAoYGFuaW1hdGlvbmAgKyBgZXh0cmVtZVN0YXRgKSAgCioqQmVycnkgQm9lc3Nlbmtvb2wqKgoKXApbdG9wXSgjdG9wKQoKIyBIeWRtb2QKSHlkcm9sb2dpY2FsIG1vZGVsbGluZyB3aXRoIGBhaXJHUmAgICAKKipLYXRpZSBTbWl0aCoqCgpcClt0b3BdKCN0b3ApCgojIEVEQQpFeHBsb3JhdG9yeSBEYXRhIEFuYWx5c2lzIGluY2x1ZGluZyBmbG93IGR1cmF0aW9uIGN1cnZlIGFuZCB0cmVuZCBhbmFseXNpcyBvbiB0aW1lLXNlcmllcyAgIAoqKlNoYXVuIEhhcnJpZ2FuKioKClwKW3RvcF0oI3RvcCkKCgojIERpc2N1c3Npb24KClBsZWFzZSBnaXZlIHVzIGZlZWRiYWNrIGF0Cjxmb250IHNpemU9IjYiPltiaXQubHkvZmVlZGJhY2tSXShodHRwczovL2JpdC5seS9mZWVkYmFja1IpPC9mb250PiAKCkZvciBkaXNjdXNzaW9ucywgcGxlYXNlIHVzZSB0aGUgCltIeWRyb2xvZ3kgaW4gUiBGYWNlYm9vayBncm91cF0oaHR0cHM6Ly93d3cuZmFjZWJvb2suY29tL2dyb3Vwcy8xMTMwMjE0Nzc3MTIzOTA5LykuICAKCg==